
METHOD FOR ESTIMATING MOVEMENT BETWEEN TWO IMAGES 



The present invention concerns a method for estimating movement between two 
numerical images. 



The movement between two successive images, h and I2, is generally defined in 
the form of a movement field associated with either of the images , h and I2 and 
constituted by movement vectors, each relating to one point of the image 
concerned. The movement vector is a two-dimensional vector representative of 
10 the difference of position between the pixel of the image h and the associated 
pixel of the image I2 relating to the same physical point of the filmed scene. 

An evaluation of movement is useful in fields for processing the image requiring a 
knowledge of movements or disparities between two images. By way of 
15 examples, it is possible to cite the following spheres of applications : 

- the compression of images : the evaluation method is used to limit the amount of 
data to code an image, the images being defined in relation to one another ; 

- the compression of data in spaces of dimension greater than 2 ; 

- video coding : the movement field defined from already coded images is then 
20 used to predict the next image ; 

- medical imagery: the method for estimating the movement between two images 
is used to conduct an analysis of the movement of the heart for example ; 

- remote monitoring : the method can be used to monitor road traffic ; 

- three-dimensional reconstruction from multi-view images: the method is used to 
25 estimate disparities between various views. 

So as to obtain this movement field, a method is known on how to break down the 
image into finished elements. These finished elements, which may for example be 
triangles or quadrangles, are determined by a meshwork whose nodes 
30 correspond to the tops of the finished elements. A movement vector is calculated 
for each node of the meshwork. Then, via the bias of an interpolation function, it is 
possible to deduce from this a movement vector for each point of the image in 
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question. The movement field is thus determined by a model of finished elements 
defining the meshwork used to partition the image into finished elements and the 
interpolation function making it possible to calculate the movement vector at any 
point of the image. 

5 

The meshwork used can be regular or irregular and needs to selected as 
sufficiently dense so as to model as best as possible the movement between the 
two images without however requiring an excessive quantity of calculation or data 
to be transmitted. This choice is made once only at the start of the method and 
10 this meshwork generally remains the same throughout the estimate. 

The calculation of the movement vectors of the nodes of the meshwork can be 
carried out according to various methods. First of all there are putting into 
correspondence or "matching' methods consisting of testing a discrete set of 

15 possible values of movement vectors for each node of the meshwork and of 
retaining the best vectors according to a given criterion. A second method known 
as a transformed method consists of using the properties of the Fourier transform 
and its extensions so as to convert the movement into a phase jump in the 
transformed space. Finally, there is a third method known as a differential method 

20 for determining the movement vectors by optimising a mathematical criterion (for 
example a quadratic error between the image and its aforesaid value with the 
movement field). This method is most frequently used for estimating movement 
with modelisation by finished elements. A conventional differential method for 
optimising movement vectors is the Gauss-Newton method. The present 

25 application concerns more particularly the movement estimation method family 
using a model of finished elements and a differential method to determine the 
movement field. 

Although widely used, this type of method does however have several drawbacks. 
30 The meshwork selected at the start of the method may prove to be inappropriate 
with respect to the semantic contents of the image, certain zones of the image 
requiring a denser meshing and others a more aerated meshing. In addition, 
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under the effect of the field of the movement vectors of the nodes of the 
meshwork, the initial meshing on the start image, for example I2, is transformed 
into a new meshing on the other image, for example U, Then pathological 
situations may then occur at the level of the new meshing, such as : 
5 - reversals of finished elements : finished elements reverse and cover others, thus 
destroying the property of partitioning of the range of the image any meshwork 
needs to verify, 

- overflowing of the peripheral nodes of the moved meshwork after applying 
movement vectors beyond the range of the image I1 : certain pixels of the image I2 
10 can be associated with pixels of the image U situated outside the range of the 
image U. This is not strictly a problem, but it may be advantageous to force the 
peripheral nodes of the moved meshwork to remain inside the range of the image 
h. 

4 0War5 AtOt> Si/inmAay of THe XlO 1/^*/ 

15 One object of the invention concerns a method for estimating movement in which 

the meshwork is optimised during estimation so as to obtain at the end of the 
method a meshwork adapted to the semantic contents of the images. To this 
effect, the finished elements are refined during the movement estimate. 

20 Another object of the invention is to improve the effectiveness of the Gauss 
Newton method so as to optimise the movement vectors of the nodes of the 
meshwork. To this effect, this optimisation is carried out on several resolution 
levels of the images. 

25 Finally, another aim of the invention is to provide a method for estimating 
movement so as to avoid aforesaid pathological situations. To this effect, the 
invention provides adding during the movement vectors optimisation step 
constraints so as to avoid these situations. 

30 Also, the invention concerns a method for estimating the movement between two 
numerical images h and I2 with luminance Yi and Y2 and intended to generate for 
each point of coordinates x,y of the image I2 a movement vector d (x,y)=(dx,dy) so 
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as to form an image I2 from the image U with luminance Y2(x,y)=Yi(x-dx,y-cly) 
which is an approximation of the image I2, characterised in that it comprises the 
following steps : 

(a) defining an initial model of finished elements comprising a meshing whose 
5 nodes are points of the image I2, a movement vector at each node of said 

meshing, and an interpolation formula for calculating the value of the movement 
vector of each point of the image b from the values of the movement vectors of 
the nodes of the mesh to which it belongs, 

(b) optimising the value of the movement vectors of the model according to a 
10 differential method, 

(c) calculating a variation E between the image I2 and the image I2 for each 
finished element or mesh, 

(d) carrying out a finer meshing on a discrete fraction of the set of finished 
elements determined according to a criterion relating to the variations E and 

15 allocating a movement vector to each new meshing node, 

(e) repeating the steps (b), (c) and (d) on the model of finished elements 
obtained at the end of the preceding step (d) until a stoppage criterion is satisfied. 

According to an improved embodiment, for each numerical image h and I2 in 
20 addition a set of R images \] is defined with a level of resolution r and luminance 
Yj' with r taking the values (0,...,R-1) and i the values 1 and 2, the images 1° and 
1° corresponding to the numerical images h and I2, the steps (b) to (e) being 
conducted for each level of resolution r from the level r=R-1 to the level r=0. 

25 Finally, according to a preferred embodiment, constraints are added to the 
movement of the finished elements at the time of optimising the movement 
vectors so as to avoid the reversal of the finished elements. According to another 
embodiment, it is also possible to introduce constraints so as to avoid the flowing 
over of the meshing obtained after applying the movement vectors beyond the 

30 sphere of the image h . 

other characteristics and advantages of the invention shall appear on a reading of 
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the following detailed description with reference to the accompanying drawings on 
which : 

- figure 1 represents a diagram of a first embodiment of the movement estimation 
method of the invention ; 

5 - figure 2 shows the step (d) of the method of the invention, and 

- figure 3 shows a diagram of an improved embodiment of the movement 
estimation method of the invention. 

Reference is first made to two numerical images h and b with respective 
10 luminance Yi and Y2. The method of the invention consists of generating for each 
point P of coordinates (x, y) in the image I2 a movement vector d(x,y)=(dx,dy). 
This vector is defined as being the vector able to construct from the image h an 
image I2 with luminance Y2 (x,y)=Yi(x-dx,y-dy) which is an approximation of I2. 
The movements are thus defined from the image h towards b. 

15 

The sought-after movement field is defined by a model of finished elements. In 
the remainder of the description, the finished elements are regarded to be 
triangles without there being any limitation of the extent of the present application 
to this form of finished elements. As a result, the model of finished elements 
20 comprises a triangular meshing, movement vectors defined in the meshing nodes, 
and an interpolation formula for calculating the movement vector of the points 
inside the triangles. 

The interpolation formula used to calculate the movement field at any point of the 
25 range of the image I2 is the following : 

If the point P of coordinates (x,y) is considered in the image I2 belonging to the 
triangle e with vertices Pj, Pj and Pr with respective coordinates (Xj.yj), (Xj,yj) and 
(Xk^yk), its movement vector is equal to 

30 H(x,y)= 2]'i'i'(>^'y)*i>yi) 

l=i.j,k 

where T,® represents a basic function associated with the triangle e 
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In the case of an affine interpolation, the ^,®(x,y) represent the barycentric 
coordinates of the point P in the triangle e with vertices Pi, Pj, Pr. These functions 
are defined by the following equation : 

f^,^ (X, y ) = a, + p,x + Y,y (x, y ) g e 

L»=i.i.k 

^r(x,y) = 0 (x,y)^e 



^iVic -x,,y= +(y: -yjx + (x,, -x = )y 

namely ^:^(x. y = — ^ 

XjVK - x,yj + x.y, - x^y, + x^y^ - x^y. 

The affine functions ^f(x,y) and ^k^(x,y) deduced from the function ^i^(x,y) by 
circularly permuting the indices i,j,k. It is also possible to use more evolved 
models of finished elements, the functions \|/ then being able to be extended to 
10 polynomials with degree n=2, but the interpolation formula of the movement 
vectors then introduces first, second derivatives, etc. A miscellaneous choice of 
models of finished elements is shown in the work "Handbook of Numerical 
Analysis" , by P.G. Ciarlet and J.L Lions, Volume 2, pages 59-99, published by 
North Holland. 

15 

According to the invention, as the movement is gradually estimated, the value of 
the movement vectors of the meshing nodes known as nodal vectors is optimised 
and the meshing is locally densified when this is necessary. Advantageously, this 
optimisation shall be carried out on several resolution levels starting with a low 
20 resolution level. 

According to a first embodiment shown on figure 1, the method of the invention 
comprises five steps referenced (a) to (e). 



According to step (a), an initial model of finished elements is defined by selecting 
points of the image I2 according to an initial triangular meshing. The nodes of the 
meshing represent the vertices of triangles (finished elements) of the model. This 
meshing can be defined in any way, for example according to the needs of the 
application or prior knowledge or of the movement already calculated between 
two preceding images of the same video sequence. If no data concerning the 
initial meshing is specified, a regular staggered meshing is used. The meshes are 
then triangles. A nil value movement vector is then associated with each node of 
the meshing. The interpolation formula previously defined is also a data element 
of the initial model. 

According to step (b), the value of the movement vectors of the model are 
optimised according to a differential method, such as the Gauss-Newton method 
or its Marquardt extension. This optimisation can be free, that is without 
constraints imposed on the possible values of the nodal vectors, or with 
constraints. The nodal vectors denote the movement vectors of the nodes of the 
meshing. Optimisation with constraints is directly linked to free optimisation and 
forms the subject of an embodiment given in detail later in the text. 

The free optimisation technique used here makes use of the advantageous 
characteristics of the Gauss-Newton method (rapid convergence of the optimum) 
and the gradient method with adaptive step (global convergence towards a local 
optimum) to resolve the linear system to be followed. This technique is an iterative 
correction of the movement vectors d(x,y) making it possible to obtain from the 
start a rough approximation of the movement. The number of iterations k for this 
optimisation of the movement vectors is either specified by the user at the start of 
the method, or depends on a threshold linked to the maximum variation between 
two consecutive values of nodal vectors for two successive iterations. Below are 
details of the Marquardt extension of the Gauss Newton optimisation method. 

The expression of the corrections dD^*^ at iteration k+1 of the movement vectors 
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according to this method is given by the following linear system : 

D*^*^ = D*'-[R'' + a.l2NrVvE*'c>-H.5D*'*^ = VE*' (1) 

with : 

- D'^*^ a column matrix of 2N elements including the components dx and dy 
of the nodal vectors on iteration k+1 , N being the number of nodes of the 
meshing in the current step : 

- D*^ a column matrix of 2N elements including the components dx and dy of 
the nodal vectors on iteration k ; 

-H = [R'^ + a,l2N] 

- I2N the identity matrix with the dimension 2N ; 



'V E'^ 
. V E"" , 



column matrix of 2N elements in which 

N elements Vx.nE*^ and N elements Vy,nE^ 

n denoting a node of the meshing and taking in turn 
the values (1 . . . N) ; with 



e€supp(n) (x,y)ee 



eGSupp(n) (x,y)GO 

where DFDk (x,y) = Y2(x.y) - Yi (x-dx.y-dy) on iteration k where sup p(n) 
represents the support of the basic function 4'„^(x,y) attached to the node 
n, that is all the triangles having the node n for a vertex ; 



where 



f j^k,x)£ j^k,xy ^ 
pk.yx Rk.yy 



a square matrix with dimension 2N 



e€supp(mn) (x,y)ee 



f ai,(x-d„y-d,r 
ax 



.H'^(x.y).^„«(x.y) 
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eGSupp(mn) (x,y)€e 



rai,(x-d,.y-d ) 



ex 



ai,( x-d,.y-dy) 

ay 



.T«(x,y).^;(x,y) 



mn mn 



e€supp(mn) (x,y)€e 



ai,(x-d,>y-dy) ^ 
dy 



.^^(x,y).^,^(x,y) 



where m and n denote nodes of the meshing and take in turn the values 
(1 ...N) and where supp(nm)=supp(n) r\ supp(m). 



- a = max JV„EN^K 



where ||^n||'s a functional norm of Tn- The two most advantageous norms 
are : 



Ml = sup(,,y^,,pp(„)|^,(x, y)| =1 or 
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1 



ZK(x,y)r 



^" ^| |sup p{n)| (x.y)^pi(n)' 
|supp(n)| denotes the cardinal number of the discrete region sup p(n). 



15 



20 



At the end of this optimisation phase, there are available N nodal vectors each 
relating to one node of the meshing. 

According to a variant embodiment of the adaptive gradient, it is possible to 
consider using a decomposition technique known as an "LDL* profile" in technical 
language so as to resolve the linear system (1 ) and accelerate the treatment. This 
technique is described in the work entitled "Matrix numerical analysis applied to 
the art of the engineer" by Theodore Lascaux, Volume 1 , pp 295-299, published 
by Masson, 1986. 
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According to an advantageous characteristic of the invention, the meshing is next 
locally refined via the division of triangles when the variation between the image I2 
and the image I2 on these triangles is too high. This is why according to the step 
(c) of the method a variation E is calculated between the image I2 and the image 
5 I2 for each triangle e. The variation E is defined as follows : 

E= ;^DFD^(x,y) 

(x.y)eo 

With DFD(x,y) = Y2 (x,y) - Yi (x-dx.y-dy) 

Of course, so as to calculate this variation for each triangle, it is necessary to 
firstly have calculated the value of the movement vectors of all the points of the 
10 image I2 by means of interpolation from the nodal vectors obtained at the end of 
step (b). 

Then in accordance with step (d), the meshing is refined on a discrete fraction of 
all the triangles of the model. This fraction is determined according to a criterion 
15 relating to the variations E previously calculated in step (c). So as to carry out this 
refining, it is possible for example to classify the triangles of the model by a 
decreasing order of their variations E and subdividing the X first triangles of this 
classification into smaller triangles. X is a predetermined fraction of the number of 
finished elements in the model, for example half. 

20 

So as to locally refine the meshing, it is also possible to compare all the variations 
E calculated in step (c) with a threshold variation which depends on the size of the 
finished element in question and of subdividing into smaller finished elements the 
finished elements whose variations E are greater than the threshold variation. 

25 

The subdivision of a triangle e into four smaller triangles is shown on figure 2. The 
triangle e is defined by the three vertices Pi, P2 and P3 having for respective 
movement vectors d^, 83, 83. So as to subdivide it into four, three new nodes P4, 
P5, Pe are defined in the middle of the three sides P1P3, P1P2, P2P3 of the triangle. 
30 Allocated to each of these three new nodes is a movement vector equal to the 
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average of the movement vectors of the two tops of the side to which it belongs, 
respectively (H,h-3^)/2, ^^+^2)^2, ^2-\-^^)/2 . The adjacent triangles to the 
triangle e whose side is P1P2, P2P3 or PiPaare then subdivided into two or three. 

5 Thus a model of finished elements is obtained whose meshing has been locally 
refined. According to step (e), steps (b), (c) and (d) are repeated at the end of the 
preceding step (d). This succession of steps is then repeated until a stoppage 
criterion is satisfied. This stoppage criterion is for example a predetermined 
number of finished elements to be reached at the end of step (d). 

10 

It is also possible to stop the method when the variations E of all these finished 
elements of the model obtained at the end of the preceding step (c) are lower 
than a threshold variation. 

15 According to an improved embodiment shown on figure 3, the steps (b) to (e) are 
carried out by depending on several resolution levels of images U and I2. The aim 
of this variant is to improve and accelerate the convergence of the calculations of 
the movement vectors. In order to achieve this, first of all for each pair of 
numerical images U and I2, a set of R images l[ is defined with a level of 

20 resolution r and luminance Y' , r taking in turn the values (R-1, R-2,...0) and i the 
values (1,2) and then steps (b) to (e) are carried out for each level of resolution r 
from the level of resolution r=R-1 to the level r=0. It is to be noted that the images 
1° and I2 correspond to the numerical images h and I2. 

25 In practice, the images are obtained by the filtering of the image Ij using a linear 
low pass filter only allowing 1/2'' of the spectral band of the image in question in 
the directions x and y, that is a pulse response filter having a pass-band 
BPr=[1/2'"'\ 1/2'""^) in the space of standardised frequencies [-1/2, 1/2]. The image 
l[ is defined by the following equation : 



30 



M M 

r(x.y)= Z ZYi(x-u.y-vKh; 

u=-M v=-M 
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The filter used is for example an approximation of an ideal filter and its pulse 
response is defined as follows : 

h;^=^ with S- ^s;, -M<n<M 

S n=-M 

and s^ = 2B.sinc(27cB,n) = 2B ^'"^^^^n 

27cBrn 

1 



B = 



2r+1 

where B and M are natural integers 
M = +00 in the ideal case. 



10 As indicated previously, this optimisation on several resolution levels is able to 
improve and accelerate the convergence of the calculations of the movement 
vectors. It is to be noted that the number of resolution levels R selected may differ 
from the number of successive refinings carried out on the meshing. 

15 According to a preferred embodiment, compactness constraints are added on 
each triangle of the model so as to prevent the triangles from reversing. 

The compactness of a triangle with vertices Pi, Pj and Pk is defined by the 
following equation : 

with C(Pi,Pj,Pk) G ]0,1[ ; and S(Pi,Pj,Pk) et P(Pi,Pj,Pk) representing 
respectively the surface and perimeter of the triangle (Pi, Pj, Pk) 
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If the compactness of a triangle is prevented from tending towards zero, it is also 
prevented from reversing. This is why, so as to avoid the reversals of triangles, 
each triangle must verify the following constraint : 
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c (p, + d, . + , p, + 3; J > K X c (p, p, p J 
<z>Kxc(p,p^,pj-c(p,+5,.,p^+dp^,p, + d;,)< 

c:>g,(D)<0 ; e = triangle(Pi,Pj,P,) 

where K is a parameter fixing the authorised compactness variation and 
D is the column vector of the movement vectors of the nodes of the model. 

According to the invention, this constraint is associated with each triangle at the 
time the movement vectors are optimised. The step for optimising the movement 
vectors amounts to a system of the type : 

ming e(^) 
^g3(^)<0 Veel 

where: 

- E(0) represents the variation between the image b and the said image I2 ; 

- g© is a constraint related to the triangle e ; 

- I is the set of the triangles of the meshing. 

So as to resolve the optimisation problems under constraints, the so-called 
increased Lagrangian technique is used. This technique is described in the work 
entitled 'Theories and algorithms" by Michel Minoux, Volume 1, pp 257-260, 
published by Dunod 1983. This technique combines two optimisation techniques : 
Lagrangian optimisation and the optimisation of external penalties. 

According to this technique, resolving the preceding system amounts to resolving 
the system without constraints according to : 

min^[^E(D)4-X:G(g3(^)?i,,rjj 

Where re is a penality element. 
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Xe is a Lagrange multiplier 

G is an increased Lagrangian determined by the equation 



G(g,(^)x,.rJ = 



Xge(D)+r,g,p) si r,>0 et ge(^)^0 
X,g,p) si r, =0 et ge(^)>0 



SI g, 



,(^)<o 



5 The constraints ge have been previously linearised by the Taylor formula to the 
order 1 : 
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The optimisation method is then the following : 

- k=0 is initialised 

- X, = 0 and r=0 are placed, X^9{"'an6 r g 9?"", and m denotes the number 
of triangles of the model 

- then the minimum 5Dk+i(A,,r) is determined so that on iteration k+1, 

- H.SD*^*^ = VE' - CV 
where Y-(X,r), C* is a matrix of xm^"" 

CV forms a matrix of the linearised constraints having for coefficients the 
following values : 

fO si ge(D)<0 



53;,ge(Pi.Pj.Pk) si P, =Pi,Pi ou P, 



sinon 



0 sinon 
where e represents the triangle with the vertices Pj, Pj and Pk. 

Then X is updated by the Uzawa algorithm and r is increased. Then the 
preceding operation is repeated until all the constraints are verified before moving 
on to iteration k+2. It is to be noted that details of the Uzawa algorithm are given 
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in the work entitled "Theories and algorithms" by Michel Minoux, Volume 1, 
published by Dunod 1983. 

According to a final embodiment, it is also possible to introduce constraints so as 
5 to avoid the flowing over of the meshing obtained after applying the movement 
vectors beyond the range of the image h. This embodiment consists of forcing the 
peripheral nodes of the meshing to remain on the edges of the image after 
applying movement vectors. In order to do this, the abscissae components 5D^*^ 
for the peripheral nodes on the left and right edges of the image I2 are cancelled 
10 on each iteration k. Similarly, the ordinate components 60*^"'^ for the upper and 
lower edges of the image I2 are cancelled on each iteration k 

In the case of a coding application, associated with the meshing is a partially 
quaternary tree. On each subdivision of the meshing (step d), an additional level 

15 is added in the tree. Each level of the tree then represents a meshing level and 
each node of the tree represents a triangle of the corresponding meshing level. 
The binary train generated during the coding step is obtained via a reading of the 
tree level by level. In this binary train, the movement vectors associated with each 
node of the tree are advantageously coded differentially with respect to the 

20 movement vectors of their father node when the latter exists and are arranged 
level by level. The corresponding decoding step consists of regenerating this tree 
from the binary train received derived from the encoder. 



